      subroutine cmfd

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'sanm.h'
      include 'setls.h'
      include 'cmfd.h'
      include 'eigv.h'
      include 'nem.h'

      real(NBF) jcor1,jcorp1,df1,chrt1,wt1,jcor2,jcorp2,df2,chrt2,wt2 &
     ,         jl1,jr1,jl2,jr2

! define dtilde, coupling coefficients
      do im=1,ng
        do j=1,tnode
          do iem=1,2
            dtilf(idm(iem,iem),1,j,im)=2*beta(iem,j-1,im)*beta(iem,j,im) &      ! left surface
                                      /(beta(iem,j-1,im)+beta(iem,j,im))      
            dtilf(idm(iem,iem),2,j,im)=2*beta(iem,j,im)*beta(iem,j+1,im) &    ! right suface
                                      /(beta(iem,j,im)+beta(iem,j+1,im))
          enddo
        enddo
! boundary condition
        dtilf(:,1,1,im)=bcm_p3(1,:)
        dtilf(:,2,j-1,im)=bcm_p3(2,:)
      enddo
! coupling coefficients
      do im=1,ng
        do j=1,tnode
          do iem=1,2
            jm1=j-1
            jp1=j+1
#define cmfd
#ifdef cmfd
! current correction factors
            dhat(idm(iem,iem),1,j,im)=-(jnet(iem,j,im)+dtilf(idm(iem,iem),1,j,im) &   
                                     *(psphi(iem,j,im)-psphi(iem,jm1,im))) &
                                     /(psphi(iem,j,im)+psphi(iem,jm1,im))
            dhat(idm(iem,iem),2,j,im)=-(jnet(iem,jp1,im)+dtilf(idm(iem,iem),2,j,im) &
                                     *(psphi(iem,jp1,im)-psphi(iem,j,im))) &
                                     /(psphi(iem,jp1,im)+psphi(iem,j,im))
#else
            dhat(idm(iem,iem),1,j,im)=0
            dhat(idm(iem,iem),2,j,im)=0
#endif
!NEM P3
#define NEMP3
#ifdef NEMP3
            btm1=beta(1,jm1,im)
            bt1=beta(1,j,im)
            btp1=beta(1,jp1,im)
            btm2=beta(2,jm1,im)
            bt2=beta(2,j,im)
            btp2=beta(2,jp1,im)            
            ppm1=psphi(1,jm1,im)
            pp1=psphi(1,j,im)
            ppp1=psphi(1,jp1,im)
            ppm2=psphi(2,jm1,im)
            pp2=psphi(2,j,im)
            ppp2=psphi(2,jp1,im)
! flux correction factors
            sffdm1=(bt1*pp1+btm1*ppm1)/(bt1+btm1)
            sffdm2=(bt2*pp2+btm2*ppm2)/(bt2+btm2)
            fcf(1,j,im)=(sflux(1,j,im)-sffdm1)/(pp1+ppm1)    
            fcf(2,j,im)=(sflux(2,j,im)-sffdm2)/(pp2+ppm2)
            if (j.eq.tnode) then
              sffdm1=(btp1*ppp1+bt1*pp1)/(btp1+bt1)
              sffdm2=(btp2*ppp2+bt2*pp2)/(btp2+bt2)
              fcf(1,jp1,im)=(sflux(1,jp1,im)-sffdm1)/(pp1+ppp1)              
              fcf(2,jp1,im)=(sflux(2,jp1,im)-sffdm2)/(pp2+ppp2)              
            endif
#endif
          enddo
        enddo
        if (bcb.eq.1) then
          dhat(1,1,1,im)=(0.375*sflux(2,1,im)-jnet(1,1,im))/sflux(1,1,im)-dtilf(1,1,1,im)
          dhat(4,1,1,im)=(0.125*sflux(1,1,im)-jnet(2,1,im))/sflux(2,1,im)-dtilf(4,1,1,im)
        endif
        if (bcu.eq.1) then
          dhat(1,2,tnode,im)=(-0.375*sflux(2,tnode+1,im)-jnet(1,tnode+1,im))/sflux(1,tnode+1,im) &
                            +dtilf(1,2,tnode,im) 
          dhat(4,2,tnode,im)=(-0.125*sflux(1,tnode+1,im)-jnet(2,tnode+1,im))/sflux(2,tnode+1,im) &
                            +dtilf(4,2,tnode,im) 
        endif        
      enddo

!      print *," "

! Current correction coefficient(dhat) adjustment-under relaxation
!#define corr
#ifdef corr
      print*,"dhat correction version"
      do im=1,ng
        do j=1,tnode
          do k=1,2
            jcor1=dhat(1,k,j,im)
            jcorp1=dhato(1,k,j,im)
            df1=dtilf(1,k,j,im)
            chrt1=(jcor1-jcorp1)/df1
!
            jcor2=dhat(4,k,j,im)
            jcorp2=dhato(4,k,j,im)
            df2=dtilf(4,k,j,im)
            chrt2=(jcor2-jcorp2)/df2
            if (abs(chrt1).gt.2.0) then
              wt1=1./(abs(chrt1)-1)
              dhat(1,k,j,im)=wt1*jcor1+(1-wt1)*jcorp1
            else if (abs(chrt2).gt.1000) then
              wt2=1./(abs(chrt2)-1)
              dhat(4,k,j,im)=wt2*jcor2+(1-wt2)*jcorp2
            endif
!            if (dhat(idm(iem,iem),k,j,im).gt.10000) then
!              print*,"heavy dhat"
!            endif
            dhato(:,k,j,im)=dhat(:,k,j,im)
          enddo
        enddo
      enddo
#endif
!
      do im=1,ng
        do j=1,tnode
!         diagonal xsec elements      
          diagf(idm(1,1),j,im)=xsr(j,im)*hz(j)
	    diagf(idm(1,2),j,im)=-2*xsr(j,im)*hz(j)
	    diagf(idm(2,1),j,im)=-2*de3*xsr(j,im)*hz(j)
	    diagf(idm(2,2),j,im)=(4*xsr(j,im)+5*xst(j,im))*de3*hz(j)
          do i=1,4
            offdiag=dtilf(i,1,j,im)+dhat(i,1,j,im)+dtilf(i,2,j,im)-dhat(i,2,j,im)
            cczf(i,1,j,im)=-dtilf(i,1,j,im)+dhat(i,1,j,im)   ! left block matrix
            cczf(i,2,j,im)=-dtilf(i,2,j,im)-dhat(i,2,j,im)   ! right block matrix
            diagf(i,j,im)=diagf(i,j,im)+offdiag    ! center block matrix
          enddo
!          dprob(j,im)=abs(diagf(1,j,im)-diagf(4,j,im))
        enddo
      enddo

      call blufac

      do im=1,ng
        do j=1,tnode
          ss=0
          do m=1,ng         
            ss=ss+avgflx(1,j,m)*xssm(j,im)%from(m)
          enddo
          src(1,j,im)=reigv*xchi(j,im)*psif(j)+ss*hz(j)
          src(2,j,im)=-2*de3*src(1,j,im)
        enddo

        call blusol(im)

        do j=1,tnode
          avgflx(1,j,im)=psphi(1,j,im)-2*psphi(2,j,im)
          avgflx(2,j,im)=psphi(2,j,im)
        enddo
      enddo

      call dumflx

! Determine the incoming partial current jin
! NEM P3 (using a pseudo flux)
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          btm2=beta(2,jm1,im)
          bt2=beta(2,j,im)
          btp2=beta(2,jp1,im)            
          ppm1=psphi(1,jm1,im)
          pp1=psphi(1,j,im)
          ppp1=psphi(1,jp1,im)
          ppm2=psphi(2,jm1,im)
          pp2=psphi(2,j,im)
          ppp2=psphi(2,jp1,im)
          sfml1=(btm1*ppm1+bt1*pp1)/(btm1+bt1)
          sfmr1=(btp1*ppp1+bt1*pp1)/(btp1+bt1)
          sfml2=(btm2*ppm2+bt2*pp2)/(btm2+bt2)
          sfmr2=(btp2*ppp2+bt2*pp2)/(btp2+bt2)
!
          sfxl1=sfml1+fcf(1,j,im)*(ppm1+pp1)
          sfxr1=sfmr1+fcf(1,jp1,im)*(ppp1+pp1) 
          sfxl2=sfml2+fcf(2,j,im)*(ppm2+pp2)
          sfxr2=sfmr2+fcf(2,jp1,im)*(ppp2+pp2) 
!
          jl1=-dtilf(1,1,j,im)*(pp1-ppm1)-dhat(1,1,j,im)*(pp1+ppm1)
          jr1=-dtilf(1,2,j,im)*(ppp1-pp1)-dhat(1,2,j,im)*(ppp1+pp1)
          jl2=-dtilf(4,1,j,im)*(pp2-ppm2)-dhat(4,1,j,im)*(pp2+ppm2)
          jr2=-dtilf(4,2,j,im)*(ppp2-pp2)-dhat(4,2,j,im)*(ppp2+pp2)          
!
          ajil(1,j,im)=0.25*sfxl1+0.5*jl1-0.1875*sfxl2
          ajir(1,j,im)=0.25*sfxr1-0.5*jr1-0.1875*sfxr2
          ajil(2,j,im)=0.4375*sfxl2+0.5*jl2-0.0625*sfxl1
          ajir(2,j,im)=0.4375*sfxr2-0.5*jr2-0.0625*sfxr1       
        enddo
        if (bcb.eq.1) ajil(1,1,im)=0
        if (bcu.eq.1) ajir(1,tnode,im)=0  
      enddo

      end subroutine
